!!!! Fortran Toolbox - (12/21/2022) !!!!
!! Author: Michael Irwin
!! Subroutines written for Fortran Code

!!!! Table of Contents !!!!
!! 1.) Linspace
!! 2.) Logspace
!! 3.) Interpolation
!! 4.) Rouwenhorst
!! 5.) Ergodic
!! 6.) Matrix Transpose
!! 7.) Matrix Inversion
!! 8.) Matrix Multiplication
!! 9.) HP Filter
!! 10.) Regression Coefficients
!! 11.) Mean
!! 12.) Standard Deviation
!! 13.) Correlation Coefficient

module Fortran_TB

use BCI_Par

implicit none

contains

!!Introduction to Code
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 1.) Linspace Subroutine

subroutine linspace(x_min, x_max, nx, xgrid)

!!!! Variables and Arrays !!!!
integer, intent(in) :: nx
real(rp), intent(in) :: x_min, x_max
real(rp), dimension(nx), intent(out) :: xgrid
integer :: x
real(rp) :: omega

!!!! Linspace Function !!!!
omega = (x_max-x_min)/(dble(nx)-one)
xgrid(1) = x_min
do x = 1,(nx-2)
	xgrid(x+1) = xgrid(x) + omega
end do
xgrid(nx) = x_max

end subroutine linspace

!! 1.) Linspace Subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 2.) Logspace Subroutine

subroutine logspace(x_min, x_max, nx, xgrid)

!!!! Variables and Arrays !!!!
integer, intent(in) :: nx
real(rp), intent(in) :: x_min, x_max
real(rp), dimension(nx), intent(out) :: xgrid
integer :: x
real(rp) :: omega

!!!! Logspace Grid !!!!
xgrid(1) = log10(x_min+one)
xgrid(nx) = log10(x_max+one)
call linspace(xgrid(1), xgrid(nx), nx, xgrid)
xgrid = 10.0_rp**xgrid
xgrid = xgrid - one

end subroutine logspace

!! 2.) Logspace Subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 3.) Interpolation Subroutine

subroutine interpolation(xgrid, xval, x_low, x_high, nx, ix, omega_x)

!!!!! Variables and Arrays !!!!
integer, intent(in) :: nx, x_low, x_high
real(rp), intent(in) :: xval
real(rp), dimension(nx), intent(in) :: xgrid
integer, intent(out) :: ix
real(rp), intent(out) :: omega_x
integer :: i_mid, i_low, i_high, diff

!!!! Interpolate !!!!
diff = 100
if ( xval >= xgrid(x_high) ) then
    ix = x_high-1
    omega_x = zero
elseif ( xval <= xgrid(x_low) ) then
    ix = x_low
    omega_x = one
else

    i_low = x_low
    i_high = x_high
    do while (diff > 1)
        i_mid = floor( (dble(i_high)+dble(i_low))/2.0_rp )
        if ( xgrid(i_mid) <= xval ) then
            i_low = i_mid
        else
            i_high = i_mid
        end if 
        diff = i_high - i_low
    end do
    ix = i_low
    omega_x = (xgrid(ix+1)-xval)/(xgrid(ix+1)-xgrid(ix))
    
end if

end subroutine interpolation

!! 3.) Interpolation Subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 4.) Rouwenhorst Subroutine

subroutine rouwenhorst(rho, sigma_eta, ne, pi_e, xgrid)

!!!! Variables and Arrays !!!!
integer :: i, m, n
integer, intent(in) :: ne
real(rp) :: p, q, sigma_e, gamma, omega
real(rp), parameter :: zero=0.0_rp, one=1.0_rp, two=2.0_rp
real(rp), intent(in) :: rho, sigma_eta
real(rp), dimension(ne), intent(out) :: xgrid
real(rp), dimension(:,:), allocatable :: p_e, pi_e2
real(rp), dimension(ne,ne), intent(out) :: pi_e

!!!! Grid of Shock Values !!!!
xgrid = zero
sigma_e = sigma_eta**two
sigma_e = sigma_e/(one-rho**two)
sigma_e = sqrt(sigma_e)
gamma = sqrt(dble(ne-1)*sigma_e**two)
xgrid(1) = -gamma
omega = (gamma+gamma)/dble(ne-1)
do m = 2,ne
	xgrid(m) = xgrid(m-1)+omega 
end do

!!!! Binomial Persistence Values !!!!
p = (rho+one)/two
q = p

!!!! Creating the Transition Matrix !!!!
i = 2
allocate(p_e(i,i))
allocate(pi_e2(i,i))
pi_e2(1,1) = p
pi_e2(1,2) = one - p
pi_e2(2,1) = one - q
pi_e2(2,2) = q

do i = 3,ne
	deallocate(p_e)
	allocate(p_e(i,i))
	p_e = zero

	!!First Block
	do m = 1,(i-1)
		p_e(m,1:i-1) = p_e(m,1:i-1) + p*pi_e2(m,1:i-1)
	end do
	p_e(i,:) = p_e(i,:) + zero
	p_e(:,i) = p_e(:,i) + zero

	!!Second Block
	do m = 1,(i-1)
		p_e(m,2:i) = p_e(m,2:i) + (one-p)*pi_e2(m,1:i-1)
	end do
	p_e(i,:) = p_e(i,:) + zero
	p_e(:,1) = p_e(:,1) + zero

	!!Third Block
	do m = 1,(i-1)
		p_e(m+1,1:i-1) = p_e(m+1,1:i-1) + (one-q)*pi_e2(m,1:i-1)
	end do
	p_e(1,:) = p_e(1,:) + zero
	p_e(:,i) = p_e(:,i) + zero

	!!Fourth Block
	do m = 1,(i-1)
		p_e(m+1,2:i) = p_e(m+1,2:i) + q*pi_e2(m,1:i-1)
	end do
	p_e(1,:) = p_e(1,:) + zero
	p_e(:,1) = p_e(:,1) + zero

	!!Assigning Value to the Transition Matrix
	deallocate(pi_e2)
	allocate(pi_e2(i,i))
	pi_e2(1,:) = p_e(1,:)
	pi_e2(i,:) = p_e(i,:)
	pi_e2(2:(i-1),:) = dble(0.5)*p_e(2:(i-1),:)

end do
pi_e = pi_e2

end subroutine rouwenhorst

!! 4.) Rouwenhorst Subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 5.) Ergodic Subroutine

subroutine ergodic(nx, pi_x, pi_erg)

!!!! Variables and Arrays !!!!
integer, intent(in) :: nx
integer :: x, x2
real(rp), parameter :: precision=0.000001_rp
real(rp) :: error
real(rp), dimension(nx) :: pi_erg2
real(rp), dimension(nx), intent(out) :: pi_erg
real(rp), dimension(nx,nx), intent(in) :: pi_x

!!!! Ergodic Distribution of Shock Values !!!!
error = 100.0_rp
pi_erg(:) = one/dble(nx)
do while (error > precision)
    pi_erg2(:) = zero
	do x = 1,nx
		do x2 = 1,nx
			pi_erg2(x2) = pi_erg2(x2) + pi_x(x,x2)*pi_erg(x)
		end do
	end do
	error = maxval(abs(pi_erg-pi_erg2))
	pi_erg = pi_erg2
end do

end subroutine ergodic

!! 5.) Ergodic Subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 6.) Matrix Tranpose 

subroutine matrix_transpose(N, M, X, Xt)

!!!! Variables and Arrays !!!!
integer, intent(in) :: N, M
real(rp), dimension(N,M), intent(in) :: X
real(rp), dimension(M,N), intent(out) :: Xt
integer :: i, j

!!!! Transpose !!!!
do i = 1,N
    do j = 1,M
        Xt(j,i) = X(i,j)
    end do
end do

end subroutine matrix_transpose

!! 6.) Matrix Transpose
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 7.) Matrix Inversion 

subroutine matrix_inversion(N, A, A_inv)

!! Referencing code by Alex G. December 2009, based on Doolittle LU factorization

!!!! Variables and Arrays (input/output) !!!!
integer, intent(in) :: N
real(rp), dimension(N,N), intent(in) :: A
real(rp), dimension(N,N), intent(out) :: A_inv
integer :: i, j, k
real(rp) :: coeff
real(rp), dimension(N) :: b, d, x
real(rp), dimension(N,N) :: L, U, A_temp

!!!! Initialize Arrays !!!!
L = zero
U = zero
b = zero
A_temp = A

!!!! Forward Elimination !!!!
do k = 1,(N-1)
    do i = (k+1),N
        coeff = A_temp(i,k)/A_temp(k,k)
        L(i,k) = coeff
        do j = (k+1),N
            A_temp(i,j) = A_temp(i,j) - coeff*A_temp(k,j)
        end do
    end do
end do

!!!! Prepare L and U Matrices !!!!
do i = 1,N
    L(i,i) = one
end do

do j = 1,N
    do i = 1,j
        U(i,j) = A_temp(i,j)
    end do
end do

!!!! Compute Columns of the Inverse Matrix !!!!
do k = 1,N
    b(k) = one
    d(1) = b(1)
    do i = 2,N
        d(i) = b(i)
        do j = 1,i-1
            d(i) = d(i) - L(i,j)*d(j)
        end do
    end do
    x(n) = d(n)/U(n,n)
    do i = (n-1),1,-1
        x(i) = d(i)
        do j = n,(i+1),-1
            x(i) = x(i) - U(i,j)*x(j)
        end do
        x(i) = x(i)/U(i,i)
    end do
    do i = 1,N
        A_inv(i,k) = x(i)
    end do
    b(k) = zero
end do

end subroutine matrix_inversion

!! 7.) Matrix Inversion
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 8.) Matrix Multiplication 

subroutine matrix_multiplication(N, Na, Nb, A, B, X)

!!!! Variables and Arrays (input/output) !!!!
integer, intent(in) :: N, Na, Nb
real(rp), dimension(Na,N), intent(in) :: A
real(rp), dimension(N,Nb), intent(in) :: B
real(rp), dimension(Na,Nb), intent(out) :: X
integer :: i, j

!!!! Matrix Multiplication !!!!
do i = 1,Na
    do j = 1,Nb
        X(i,j) = dot_product(A(i,:),B(:,j))
    end do
end do

end subroutine matrix_multiplication

!! 8.) Matrix Multiplication
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 9.) HP Filter

subroutine HP_Filter(N, lambda, x_vec, x_tilde)

!!!! Variables and Arrays (input/output) !!!!
integer, intent(in) :: N
real(rp), intent(in) :: lambda
real(rp), dimension(N), intent(in) :: x_vec
real(rp), dimension(N), intent(out) :: x_tilde
integer :: t
real(rp), dimension(N) :: x_hat
real(rp), dimension(N,N) :: A, A_inv

!!!! Row 1 !!!!
A(:,:) = zero
A_inv(:,:) = zero
A(1,1) = one + lambda
A(1,2) = -2.0_rp*lambda
A(1,3) = lambda

!!!! Row 2 !!!!
A(2,1) = -2.0_rp*lambda
A(2,2) = one + 5.0_rp*lambda
A(2,3) = -4.0_rp*lambda
A(2,4) = lambda

!!!! Rows 3 through N-2 !!!!
do t = 3,N-2
    A(t,t-2) = lambda
    A(t,t-1) = -4.0_rp*lambda
    A(t,t) = one + 6.0_rp*lambda
    A(t,t+1) = -4.0_rp*lambda
    A(t,t+2) = lambda
end do

!!!! Row N-1 !!!!
A(N-1,N-3) = lambda
A(N-1,N-2) = -4.0_rp*lambda
A(N-1,N-1) = one + 5.0_rp*lambda
A(N-1,N) = -2.0_rp*lambda

!!!! Row N !!!!
A(N,N-2) = lambda
A(N,N-1) = -2.0_rp*lambda
A(N,N) = one + lambda

!!!! Create Inverted Matrix !!!!
call matrix_inversion(N, A, A_inv)

!!!! Create Trend !!!!
call matrix_multiplication(N, N, 1, A_inv, x_vec, x_hat)

!!!! Detrend Data !!!!
x_tilde(:) = x_vec(:) - x_hat(:)


end subroutine HP_Filter

!! 9.) HP Filter 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 10.) Regression Coefficients

subroutine regression_coefficients(N, Nx, y_vec, X, b_vec)

!! Solves: (X'X)^{-1} X'y = b

!!!! Input/Output Variables !!!!
integer, intent(in) :: N, Nx
real(rp), dimension(N), intent(in) :: y_vec
real(rp), dimension(N,Nx), intent(in) :: X
real(rp), dimension(Nx), intent(out) :: b_vec
real(rp), dimension(Nx,N) :: XT, XTXXT
real(rp), dimension(Nx,Nx) :: XTX, XTX_inv

!!!! Calculate Regression Coefficients !!!!
call matrix_transpose(N, Nx, X, XT)
call matrix_multiplication(N, Nx, Nx, XT, X, XTX)
call matrix_inversion(Nx, XTX, XTX_inv)
call matrix_multiplication(Nx, Nx, N, XTX_inv, XT, XTXXT)
call matrix_multiplication(N, Nx, 1, XTXXT, y_vec, b_vec)

end subroutine regression_coefficients

!! 10.) Regression Coefficients
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 11.) Mean

subroutine average(N, vector, mean)

!!!! Variables and Arrays !!!!
integer, intent(in) :: N
real(rp), dimension(N), intent(in) :: vector
real(rp), intent(out) :: mean

!!!! Calculate Mean !!!!
mean = sum(vector(1:N))/dble(N)

end subroutine average

!! 11.) Mean
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 12.) Standard Deviation

subroutine standard_deviation(N, mean, vector, sd)

!!!! Variables and Arrays !!!!
integer, intent(in) :: N
real(rp), intent(in) :: mean
real(rp), dimension(N), intent(in) :: vector
real(rp), intent(out) :: sd
integer :: i

!!!! Calculate Standard Deviation !!!!
sd = zero
do i = 1,N
    sd = sd + ( (vector(i) - mean)**2.0_rp )/dble(N)
end do
sd = sqrt(sd)

end subroutine standard_deviation

!! 12.) Standard Deviation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 13.) Correlation Coefficient

subroutine correlation_coefficient(N, mean_x, mean_y, x_vec, y_vec, corr)

!!!! Variables and Arrays !!!!
integer, intent(in) :: N
real(rp), intent(in) :: mean_x, mean_y
real(rp), dimension(N), intent(in) :: x_vec, y_vec
real(rp), intent(out) :: corr
integer :: i
real(rp) :: numerator, denominator, x_diff, y_diff

!!!! Calculate Correlation Coefficient !!!!
numerator = zero
denominator = zero
do i = 1,N
    numerator = numerator + ( x_vec(i) - mean_x )*( y_vec(i) - mean_y )
    x_diff = x_diff + (x_vec(i)-mean_x)**2.0_rp
    y_diff = y_diff + (y_vec(i)-mean_y)**2.0_rp
end do
denominator = sqrt(x_diff*y_diff)
corr = numerator/denominator

end subroutine correlation_coefficient

!! 13.) Correlation Coefficient
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

end module Fortran_TB